Figure 1. Image of Siberian hamster Phodopus sungorus. Image attribution: Salzgeber, P. (2007). Phodopus sungorus. In Wikimedia Commons. https://commons.wikimedia.org/wiki/File:PhodopusSungorus_1.jpg
Question: How does early development (i.e., maternal
environment) affect whether offspring vary in their social behaviors
during the course of a long-bout social stressor?
Hypothesis: Maternal stress during offspring
development has sex-specific effects on consistency in individual
behavior during a long-bout social stressor
Predictions:
1. Female offspring from Stress Only
mothers will display increased aggression throughout a long-bout
stressor.
2. Other offspring will not display consistent behaviors
throughout the long-bout stressor.
Figure 2. Maternal treatment graphic depicting the four treatment groups.
Figure 3.Offspring in late adolescence (51-56 days old) were exposed to a resident-intruder trial as above. The intruder was placed in the focal individual’s cage for 15 minutes. Time scored is shown by red line.
Figure 4. Actual screenshots of video data depicting the four social behaviors of interest.
Below is the code used to clean the and prepare the data for analyses. Several columns were removed from the data frame because they were not used in this project (these were columns imported into the excel sheet from the original project in 2019).
# Load data, clean names, view structure ####
hammies <- read_csv("../final_project/IU106_Social Behavior First and Last Five Minutes_Formatted_6-6-23.csv")
hammies <- janitor::clean_names(hammies)
names(hammies)
str(hammies)
# Remove columns that were not used in this project ####
columns_removed <- c("id", "aggression_number_assigned", "comp1", "comp2", "comp3",
"aggression_score_freq", "aggression_score_duration")
hammies <- dplyr::select(hammies, -columns_removed)
# Create trial timing column, set appropriate columns as factors ####
tidy_hammies <- hammies %>%
mutate(trial_timing = as.factor(first_or_last_five_minutes)) %>%
mutate(treatment = as.factor(maternal_treatment)) %>%
mutate(juvenile_id = as.factor(juvenile_id)) %>%
dplyr::select(-first_or_last_five_minutes)
# Fix spelling of all "occurence" columns ####
tidy_hammies <- tidy_hammies %>%
mutate(attack_number_of_occurrences = attack_number_of_occurences) %>%
mutate(chase_number_of_occurrences = chase_number_of_occurences) %>%
mutate(groom_number_of_occurrences = groom_number_of_occurences) %>%
mutate(investigation_number_of_occurrences = investigation_number_of_occurences) %>%
mutate(jump_number_of_occurrences = jump_number_of_occurences) %>%
dplyr::select(-contains("occurences"))
str(tidy_hammies)
We analyzed the number (i.e., frequency) and total duration of
behavioral bouts in the first and last five minutes (t0-t5 vs. t10-t15).
For each model, we first tested for a three-way interaction between
maternal treatment, offspring sex, and trial timing. If no significant
three-way interaction was detected, we tested for a two-way interaction
between offspring sex and maternal treatment. Only significant
interactions were retained in the model; non-significant interaction
terms were removed.
All models had two random effects: intruder identity and resident (offspring) identity, to control for individual differences in behavior. All frequency models (GLMM’s) used the bobyqa optimizer from the lme4 package for model fitting.
| Response Variable | Distribution | Fixed Effects |
|---|---|---|
| Attack Frequency | Poisson, Log link function | Offspring sex, maternal treatment, trial timing, SI- CORT, offspring weight, offspring sex*maternal treatment interaction |
| Attack Duration | Gaussian, Identity link function | Offspring sex, maternal treatment, trial timing, SI-CORT, offspring weight |
| Investigation Frequency | Poisson, Log link function | Offspring sex, maternal treatment, trial timing, SI-CORT, offspring weight, offspring sexmaternal treatmenttrial timing |
| Investigation Duration | Gaussian, Identity link function | Offspring sex, maternal treatment, trial timing, SI-CORT, offspring weight |
| Escape Frequency | Poisson, Log link function | Offspring sex, maternal treatment, trial timing, SI-CORT, offspring weight, offspring sexmaternal treatmenttrial timing |
| Escape Duration | Gaussian, Identity link function | Offspring sex, maternal treatment, trial timing, SI-CORT, offspring weight |
| Grooming Frequency | Poisson, Log link function | Offspring sex, maternal treatment, trial timing, SI-CORT, offspring weight |
| Grooming Duration | Gaussian, Identity link function | Offspring sex, maternal treatment, trial timing, SI-CORT, offspring weight, offspring sex* maternal treatment |
Our first model analyzes the effects of maternal treatment, offspring sex, and video timing on the frequency of aggressive behaviors. There was no significant three-way interaction, but we did detect a significant interaction between maternal treatment and offspring sex. Below is the code for model A.1.
modelA.1 <- tidy_hammies %>%
glmer(attack_number_of_occurrences ~ sex*treatment+trial_timing+scale(si_cort)
+weight_at_trial+(1|intruder_id)+(1|juvenile_id),family=poisson,
glmerControl(optimizer="bobyqa"),data=.)
broom.mixed::tidy(modelA.1) %>%
kable("pipe", digits = 3, caption = "**Table 2.** Summary of model A.1, analyzing the interactive effects of maternal treatment and offspring sex on frequency of attack behaviors in the first and last five minutes of the resident-intruder trial.") %>%
kable_styling() %>%
kableExtra::row_spec(.,c(4:6, 10), bold=TRUE) %>%
kableExtra::column_spec(column = 1:7, width = "200px")
| effect | group | term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|---|---|
| fixed | NA | (Intercept) | 0.589 | 0.706 | 0.834 | 0.404 |
| fixed | NA | sexMale | 0.430 | 0.373 | 1.153 | 0.249 |
| fixed | NA | treatmentAntibiotic Only | 0.604 | 0.388 | 1.555 | 0.120 |
| fixed | NA | treatmentControl | -2.299 | 1.050 | -2.190 | 0.029 |
| fixed | NA | treatmentStress Only | 0.852 | 0.384 | 2.221 | 0.026 |
| fixed | NA | trial_timingLast | -0.321 | 0.153 | -2.103 | 0.035 |
| fixed | NA | scale(si_cort) | 0.023 | 0.101 | 0.228 | 0.819 |
| fixed | NA | weight_at_trial | 0.014 | 0.022 | 0.645 | 0.519 |
| fixed | NA | sexMale:treatmentAntibiotic Only | -0.551 | 0.485 | -1.137 | 0.256 |
| fixed | NA | sexMale:treatmentControl | 2.323 | 1.106 | 2.101 | 0.036 |
| fixed | NA | sexMale:treatmentStress Only | -0.784 | 0.482 | -1.626 | 0.104 |
| ran_pars | juvenile_id | sd__(Intercept) | 0.145 | NA | NA | NA |
| ran_pars | intruder_id | sd__(Intercept) | 0.000 | NA | NA | NA |
Model A.1 is a poisson mixed model fitted to predict number of attack
bouts.
tidy_hammies %>%
ggplot(aes(x=interaction(trial_timing, sex),
y = attack_number_of_occurrences))+
geom_boxplot(aes(fill=treatment))+
labs(x= "Trial timing",
y = "Number of attacks",
fill = "Maternal treatment")+
theme(axis.line = element_line(color="black"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank())
Figure 5. Plot of aggression frequency by sex and trial timing.
We then used the lsmeans function from the emmeans R package to conduct pairwise comparisons of the maternal treatment and sex interactions. lsmeans uses a corrected p-value since to minimize the likelihood of a type I error, since we are performing a large number of comparisons. After running the lsmeans, none of our pairwise comparisons are statistically significant using a p-value of 0.05.
#pairwise comparisons with corrected p-value
lsmeans(modelA.1, pairwise~treatment, type = "response") # because there was a significant 2-way interaction, we are more interested in the interaction below.
lsmeans(modelA.1, pairwise~treatment*sex, type = "response")
lsmeans_output <- lsmeans(modelA.1, pairwise~treatment*sex, type = "response")
lsmeans_contrasts <- lsmeans_output$contrasts
lsmeans_contrasts_df <- as.data.frame(lsmeans_contrasts)
lsmeans_contrasts_df %>%
kable("pipe", digits = 3, caption = "**Table 3.** lsmeans comparisons for attack frequency.") %>%
kable_styling() %>%
kableExtra::row_spec(.,14, bold=TRUE) %>%
kableExtra::column_spec(column = 1:7, width = "200px")
| contrast | estimate | SE | df | null | z.ratio | p.value |
|---|---|---|---|---|---|---|
| (Antibiotic + Stress Female) / Antibiotic Only Female | 0.547 | 0.212 | Inf | 1 | -1.555 | 0.777 |
| (Antibiotic + Stress Female) / Control Female | 9.964 | 10.463 | Inf | 1 | 2.190 | 0.358 |
| (Antibiotic + Stress Female) / Stress Only Female | 0.427 | 0.164 | Inf | 1 | -2.221 | 0.339 |
| (Antibiotic + Stress Female) / (Antibiotic + Stress Male) | 0.651 | 0.243 | Inf | 1 | -1.153 | 0.945 |
| (Antibiotic + Stress Female) / Antibiotic Only Male | 0.617 | 0.217 | Inf | 1 | -1.370 | 0.871 |
| (Antibiotic + Stress Female) / Control Male | 0.635 | 0.219 | Inf | 1 | -1.316 | 0.893 |
| (Antibiotic + Stress Female) / Stress Only Male | 0.608 | 0.231 | Inf | 1 | -1.307 | 0.897 |
| Antibiotic Only Female / Control Female | 18.225 | 19.168 | Inf | 1 | 2.760 | 0.105 |
| Antibiotic Only Female / Stress Only Female | 0.780 | 0.273 | Inf | 1 | -0.710 | 0.997 |
| Antibiotic Only Female / (Antibiotic + Stress Male) | 1.190 | 0.445 | Inf | 1 | 0.465 | 1.000 |
| Antibiotic Only Female / Antibiotic Only Male | 1.129 | 0.389 | Inf | 1 | 0.352 | 1.000 |
| Antibiotic Only Female / Control Male | 1.161 | 0.395 | Inf | 1 | 0.440 | 1.000 |
| Antibiotic Only Female / Stress Only Male | 1.112 | 0.418 | Inf | 1 | 0.282 | 1.000 |
| Control Female / Stress Only Female | 0.043 | 0.045 | Inf | 1 | -2.975 | 0.059 |
| Control Female / (Antibiotic + Stress Male) | 0.065 | 0.067 | Inf | 1 | -2.648 | 0.139 |
| Control Female / Antibiotic Only Male | 0.062 | 0.064 | Inf | 1 | -2.709 | 0.120 |
| Control Female / Control Male | 0.064 | 0.066 | Inf | 1 | -2.672 | 0.131 |
| Control Female / Stress Only Male | 0.061 | 0.063 | Inf | 1 | -2.709 | 0.120 |
| Stress Only Female / (Antibiotic + Stress Male) | 1.525 | 0.588 | Inf | 1 | 1.095 | 0.958 |
| Stress Only Female / Antibiotic Only Male | 1.447 | 0.485 | Inf | 1 | 1.103 | 0.956 |
| Stress Only Female / Control Male | 1.489 | 0.513 | Inf | 1 | 1.155 | 0.944 |
| Stress Only Female / Stress Only Male | 1.425 | 0.529 | Inf | 1 | 0.955 | 0.980 |
| (Antibiotic + Stress Male) / Antibiotic Only Male | 0.949 | 0.285 | Inf | 1 | -0.176 | 1.000 |
| (Antibiotic + Stress Male) / Control Male | 0.976 | 0.308 | Inf | 1 | -0.077 | 1.000 |
| (Antibiotic + Stress Male) / Stress Only Male | 0.935 | 0.299 | Inf | 1 | -0.212 | 1.000 |
| Antibiotic Only Male / Control Male | 1.029 | 0.307 | Inf | 1 | 0.096 | 1.000 |
| Antibiotic Only Male / Stress Only Male | 0.985 | 0.284 | Inf | 1 | -0.052 | 1.000 |
| Control Male / Stress Only Male | 0.958 | 0.317 | Inf | 1 | -0.131 | 1.000 |
Below is the code my PI initially used to check model assumptions and check for overdispersion.
# * * Checking model assumptions ------------------------------------------
qqnorm(residuals(modelA.1))
qqline(residuals(modelA.1))
shapiro.test(residuals(modelA.1))
# * * Check overdispersion ------------------------------------------------
overdisp_fun <- function(modelA.1) {
rdf <- df.residual(modelA.1)
rp <- residuals(modelA.1,type="pearson")
Pearson.chisq <- sum(rp^2)
prat <- Pearson.chisq/rdf
pval <- pchisq(Pearson.chisq, df=rdf, lower.tail=FALSE)
c(chisq=Pearson.chisq,ratio=prat,rdf=rdf,p=pval)
}
overdisp_fun
overdisp_fun(modelA.1)
sum(residuals(modelA.1,type="pearson")^2)/(nrow(tidy_hammies)-length(fixef(modelA.1)-1))
vif(modelA.1)
Using the performance package we learned about in class, I simplified the overdispersion code, and used check_model to visualize vif and overdispersion.
check_model(modelA.1)
Figure 6. Output of performance::check_model function for model A.1.
performance::check_overdispersion(modelA.1)
## # Overdispersion test
##
## dispersion ratio = 1.341
## Pearson's Chi-Squared = 46.932
## p-value = 0.086
We fitted a linear mixed model for the duration of attack behaviors in the first and last five minutes of the resident-intruder trial. The data were square-root transformed to assume a normal distribution. There were no significant interactions to predict attack behavior duration.
modelA.2 <- lmer(sqrt(attack_total_duration)~ sex+treatment+trial_timing+
scale(si_cort)+weight_at_trial+(1|intruder_id)+(1|juvenile_id),
data = tidy_hammies)
#summary(modelA.2)
broom.mixed::tidy(modelA.2) %>%
kable("pipe", digits = 3) %>%
kable_styling() %>%
kableExtra::column_spec(column = 1:7, width = "200px")
| effect | group | term | estimate | std.error | statistic | df | p.value |
|---|---|---|---|---|---|---|---|
| fixed | NA | (Intercept) | 3.608 | 1.390 | 2.597 | 17.398 | 0.019 |
| fixed | NA | sexMale | 0.749 | 0.452 | 1.656 | 17.000 | 0.116 |
| fixed | NA | treatmentAntibiotic Only | 0.324 | 0.600 | 0.540 | 17.000 | 0.596 |
| fixed | NA | treatmentControl | -0.256 | 0.586 | -0.437 | 17.000 | 0.668 |
| fixed | NA | treatmentStress Only | 0.397 | 0.631 | 0.630 | 17.000 | 0.537 |
| fixed | NA | trial_timingLast | -0.009 | 0.299 | -0.031 | 23.000 | 0.976 |
| fixed | NA | scale(si_cort) | 0.180 | 0.251 | 0.717 | 17.000 | 0.483 |
| fixed | NA | weight_at_trial | -0.069 | 0.043 | -1.584 | 17.000 | 0.132 |
| ran_pars | juvenile_id | sd__(Intercept) | 0.677 | NA | NA | NA | NA |
| ran_pars | intruder_id | sd__(Intercept) | 0.000 | NA | NA | NA | NA |
| ran_pars | Residual | sd__Observation | 1.034 | NA | NA | NA | NA |